IMPLICIT-EXPLICIT RUNGE-KUTTA SCHEMES FOR HYPERBOLIC SYSTEMS 
AND KINETIC EQUATIONS IN THE DIFFUSION LIMIT 
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Abstract. We consider Implicit- Explicit (IMEX) Rungc-Kutta (R-K) schemes for hyperbolic systems with stiff relax- 
ation in the so-called diffusion limit. In such regime the system relaxes towards a convection-diffusion equation. The first 
objective of the paper is to show that traditional partitioned IMEX R-K schemes will relax to an explicit scheme for the limit 
equation with no need of modification of the original system. Of course the explicit scheme obtained in the limit suffers from 
the classical parabolic stability restriction on the time step. The main goal of the paper is to present an approach, based 
on IMEX R-K schemes, that in the diffusion limit relaxes to an IMEX R-K scheme for the convection-diffusion equation, 
in which the diffusion is treated implicitly. This is achieved by an original reformulation of the problem, and subsequent 
application of IMEX R-K schemes to it. An analysis on such schemes to the reformulated problem shows that the schemes 
reduce to IMEX R-K schemes for the limit equation, under the same conditions derived for hyperbolic relaxation [8]. Several 
numerical examples including neutron transport equations confirm the theoretical analysis. 
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1. Introduction. The development of numerical methods to solve hyperbolic systems in diffusive 
regimes has been a very active area of research in the last years (see for example pSl [5T1 134j ). 

Classical fields of applications involve diffusion in neutron transport HSl [111 HSl 132] , drift-diffusion 
limit in semiconductors }271 130] and incompressible Navier-Stokes limits in rarefied gas dynamic f29j . A 
strictly related field of research concerns the construction of schemes for the compressible Navier-Stokes 
limit (see [3] and the references therein). In such physical problems, the scaling parameter (mean free 
path) may differ in several orders of magnitude from the rarefied regimes to the diffusive regimes, and it 
is desirable to develop a class of robust numerical schemes that can work uniformly with respect to this 
parameter. 

A prototype hyperbolic system of conservation laws with diffusive relaxation that we will use to 
illustrate the subsequent theory is the following, [25l [27l [34] 

ut + v^ = 0, 

1 1 (1-1) 

+ -2Piu)x = -—{v- q{u)) , 

where p'{u) > 0. System is hyperbolic with two distinct real characteristics speeds ±y^p'{u)/e. 

In the small relaxation limit, e 0, the behavior of the solution to is, at least formally, governed 
by the convection-diffusion equation 

Ut + q{u)x = p{u)xx, 

(1.2) 

V = q{u) -p{u)x- 

The so called subcharacteristic condition [T5] for system (|l.ip becomes 

k'(«)P<^, (1.3) 

and is naturally satisfied in the limit £ — > 0. 

About the boundary conditions for system (jl.ip . we have to specify the domain in which we solve 
the problem. In a finite domain x G [a, &], one can use periodic boundary conditions, or can assign one 
independent condition at each boundary, since the two characteristic velocities have opposite sign. 

In practice, we shall assign two conditions at each boundary, one independent and the other compat- 
ible with the equations. Furthermore, we shall choose boundary conditions which are compatible to the 
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limit solution as e — 0. For example for system if we set v = q{u) at x ~ a and x = b, compatibility 

with system (jl.ip requires 

e^q' {u)vx = p' {u)ux, X = a, b. 

Such condition becomes = in the limit e — > 0. In the sections on the numerical tests we shall specify 
the boundary conditions we use in each test. 

In general, numerical approaches that work for hyperbolic system with stiff relaxation terms do not 
apply directly in the diffusive scaling since in these systems we have the presence of multiple time-scales. 

In fact, together with the stiff relaxation term we have a stiff convection term that contributes to the 
asymptotic diffusive behavior. Then special care must be taken to ensure that the schemes possess the 
correct zero-relaxation limit, in the sense that the numerical scheme applied to system should be a 
consistent and stable scheme for the limit system ()1.2|) as the parameter s approaches zero independently 
of the discretization parameters. A notion usually referred to as asymptotic preservation. For a nice 
survey on asymptotic preserving scheme for various kinds of systems see, for example, the review paper 
by Shi Jin ([23]). Furthermore, a different approach to the derivation of asymptotic preserving schemes 
is described in the review by Pierre Degond |17j . In the case of Boltzmann kinetic equations we also refer 
to the recent review by two of the authors [37]. 

IMEX Runge-Kutta (R-K) schemes [U [5l IH [HI [S^ represent a powerful tool for the time discretization 
of such stiff systems. Unfortunately, since the characteristic speed of the hyperbolic part is of order 1/e, 
standard IMEX R-K schemes developed for hyperbolic systems with stiff relaxation [5] become useless 
in such parabolic scaling, because the CFL condition would require At — 0{eAx). Of course, in the 
diffusive regime where s < Ax, this is very restrictive since for an explicit method a parabolic condition 
At C(Aa;2) should suffice. 

Most previous works on asymptotic preserving schemes for hyperbolic systems and kinetic equations 
with diffusive relaxation focus on schemes which in the limit of infinite stiffness become consistent explicit 
schemes for the diffusive limit equation |25| [211 [3TJ [32l [34] . Such schemes have been derived by splitting 
the stiff hyperbolic part into an explicit (non-stiff) term, and an implicit (stiff) term. Here we show that 
by applying partitioned IMEX R-K schemes, in which the stiffness is associated with the variable and not 
with the operator, one obtains IMEX R-K schemes that naturally relax to the explicit scheme applied to 
the limit convection-diffusion equation. All these explicit schemes clearly suffer from the usual stability 
restriction At = 0{Ax^). 

In this paper we present a general methodology to overcome such stability restriction which applies 
to a broad class of problems. The idea is to reformulate problem by properly combining the limiting 
diffusion flux with the convective flux. This allows to construct a class of IMEX R-K schemes that 
work with high order accuracy in time and that, in the diffusion limit (i.e. when e — > 0), originate an 
IMEX method for the limiting convection-diffusion equation (|1.2p . Other reformulations whose goal is 
to obtain asymptotic-preserving methods have been proposed in |251 124j . Schemes that avoid such time 
step restriction and provide fully implicit solvers in the case of transport equations have been analyzed 
in [g. 

Our new approach allows a hyperbolic CFL condition At = 0{Ax) independent of e when applied to 
(jl.ip in all regimes. The aim of this paper is to derive and analyze different types of IMEX R-K schemes 
when applied to the reformulated problem in the stiff regime (e — )■ 0). 

The rest of the paper is organized as follows. The next section is devoted to partitioned IMEX R-K 
schemes. It is shown that they relax to the explicit scheme applied to the convection diffusion limit. In 
Section [3] the new approach is introduced and analyzed. In particular, following [8], we prove that under 
suitable assumptions the IMEX R-K schemes are consistent with the diffusion limit. The analysis is based 
on a power series expansion in e of the exact and numerical solution. It is shown that, at lowest order 
in e, the model system is a set of differential-algebraic equations of index 1, i.e. it can be transformed 
into a set of ordinary differential equations by one time differentiation. Compatibility between exact and 
numerical solution at different orders in e introduces additional order conditions on the coefficients of the 
IMEX schemes. To the lowest order such conditions are referred as index 1 order conditions. 

After a short section on space discretization obtained by conservative finite difference schemes, in 
Section [5] we report several numerical examples and tests. In Section [6] we consider one-dimensional 
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neutron transport equation and present several numerical results and comparison with schemes available 
in the literature. Additional technical material is given in the separate Appendices. 



2. Partitioned IMEX R-K schemes. The first observation is that in system (|l.ip the stiffness is 
naturally associated to the variable v rather then to some operator. The system has the structure of a 
singular perturbation problem |21| , and it can be treated by a partitioned R-K scheme in which the first 
equation is treated explicitly and the second implicitly 

ut = -Vx, (Exphcit) 

(2.1) 

e'^vt = -{p{u)x + v ~ q{u)). (Implicit) 

This approach has been used, for example, in ^51 [55] . 

By using a method of lines approach (MOL) , we discretize system p.ip in space by a uniform mesh 
{^i}i=i ^^'^ Ui{t) sa u{xi,t), Vi{t) w v{xi,t). We obtain a large system of ODE's 

Ut = TiV), 

(2.2) 

e^Vt = G{U) - V, 

with U{t) = (C/i(<),C/2(t),...,t/jv(t))^ e and V{t) ^ (Fi(i), ^2(0, Viv(t))^ G R^, where F{V) ^ 
-DV and G{U) = Q{U) - Dp{U). Here DV and Dp{U) (with a shght abuse of notation) denote the 
discretization of the convective terms v^, p{u)x, while Q{U) represents the discretization of the term q{u). 

As we shall see, if the implicit scheme is L-stable [21], in the limit e ^ the IMEX R-K scheme will 
relaxe to the explicit scheme applied to the limit equation 

Ut HU), (2.3) 

where T{U) = T{G{U)). 

For example, implicit-explicit Euler scheme applied to system p.2p gives 

Un+l = Un + At^(K) 

e^Vn + MG{Un+l) 



Vr 



n+1 



At 



where we have discretized the interval of integration by a time mesh {inj^J^i and C7„ « U{tn). As e — ?> 0, 
Vn+i = G(C/„+i) and therefore Un+i = C/„ + AtJ'(C/„). 

In the case p{u) = u and q{u) = 0, the method would relax to explicit Euler scheme applied to the 
diffusion equation, thus suffering the usual parabolic CFL stability condition At ^ Ax^/2. This approach 
will be denoted as partitioned approach. 

2.1. Classification of IMEX R-K schemes. IMEX R-K schemes have been widely used in the 
literature to treat problems that contain both stiff and non stiff terms [TJ [TT] [3S]. The stiff terms 
are treated implicitly, while the non stiff terms are treated explicitly, thus lowering the computational 
complexity of the scheme. 

Usually such a scheme is characterized by the s x s matrices A = (flij), A = (oy) and the vectors 6, 
G R'*, and can be represented by a double tableau in the usual Butcher notation 



c 


A 


C 


A 






6^ 



The coefficients c and c are used if the right hand side depends explicitly on time. We assume that they 
satisfy the usual relation 



^fiy, Q^^Oij. (2.4) 



Matrix A is lower triangular with zero diagonal, while matrix A is lower triangular, i.e. the implicit 
scheme is a diagonally implicit Runge-Kutta (DIRK). This choice guarantees that the term T{V) in (|2.2p 
is always explicitly evaluated. 

IMEX R-K schemes presented in the literature can be classified in two main different types charac- 
terized by the structure of the matrix A = (aij)f of the implicit scheme. 

Definition 2.1. We call an IMEX R-K method of type A (see f36^ } if the matrix A e R"^'' is 
invertible. 

Definition 2.2. We call an IMEX R-K method of type CK (see JTTj) if the matrix A e ^ can 
be written as 




with a G K.^""-'^) and the submatrix A € M.^'^^^^ ^ (s-i) invertible. In the special case a = the scheme is 
said to be of type ARS (see fl^). 

We note that schemes CK are very attractive because they allow some simplifying assumptions, that 
make order conditions easier to treat, therefore permitting the construction of higher order IMEX R-K 
schemes. On the other hand, schemes of type A are more amenable to a theoretical analysis, since the 
matrix A of the implicit scheme is invertible. This is why we start our analysis with the latter schemes. 

2.2. Analysis of IMEX schemes for the partitioned approach. Now as an example we perform 
the analysis of type A scheme when applied to system (|2.ip . A similar analysis is possible also for CK 
schemes. We will restrict our analysis to the limit case e — 0. 

Applying an IMEX R-K scheme to system ()2.2p we obtain 

s 

Un+i = Un + AtJ2hJ'iVk), 

k=l 

s 

e^Vn+i = eV„ + At ^ 6fc(G(Wfc) 

k=l 

for the numerical solution and 

fe-i 

Uk = Un+AtY,^kjJ'{Vj) 

i=i 

k 

e^Vk = eV„ + AtJ2akj{g{Kj) 

for the internal stages. 

By Definition 12.11 and A invertible we obtain from the second equation in (|2.5p 

k 

At{G{Uk) - Vk) = ^ ^fc,(V, - Vn). (2.7) 
i=i 

From now on, iOkj are the elements of the inverse matrix A^^ . Now inserting (j2.7p into the numerical 
solution Vn+i and setting e = 0, we get 

s 

Vn+l = n{^)Vn +AtJ2 hiOkjVj. (2.8) 

fe=l 

Here we denoted by 



(2.5) 



Vfc), 



(2.6) 



7^(cx)) = 1 - b^A-H = lim 7^(z), 



where TZ{z) is the stability function of the imphcit scheme defined by (see |21| . Sect. IV. 3) 

n{z) = l + zb^{I - zA)-^l, (2.9) 
with b'^ = (61, . . . , 6s) and 1 (I, . . . , I)^. 

By (12171) we get Vk = G{Uk) when e = 0. This yields that f{Uk) = T{G{Uk)), and we obtain 

S 

C/„+i = t^„ + At^6fc^(Z^fc) 
fe=i 

with 

j=i 

internal stages and k = 1, . . . , s. This represents the explicit scheme of the starting IMEX R-K one of 
type A applied to the limit equation (|2.3p obtained by (|2.2p when e — >• 0. As particular case, if p{u) = u 
and q{u) = 0, this is the explicit scheme applied to the limit diffusion equation under the usual parabolic 
stability restriction (At ^ Ax^/2). 

3. Overcoming parabolic stiffness. In order to overcome such stability restriction, we reformulate 
system p.ip as the equivalent system 

ut = -{v + iip{u)x)x + tJ-p{u)xx, 

(3.1) 

e^vt = -p{u)x - w + q{u), 

where the term ^p{u)xx has been added and subtracted to the first equation in (jl.ip . Here ^ = ^{e) € [0, 1] 
is a free parameter such that ^(0) = 1. The idea is that, since the quantity v + p{u)x is close to q{u) as 
e — >■ 0, the first term on the right hand side can be treated explicitly in the first equation, while the term 
p{u)xx will be treated implicitly. This can be done naturally by using an Implicit-Explicit approach, as 
we will explain later. Let us point out that the choice /i = 1, as shown in Appendix 18.11 for a first order 
implicit-explicit scheme, guarantees the largest stability region of the method. 

Next we will study the behavior of the different IMEX R-K schemes when applied to system p.ip in 
the diffusion limit. In particular we will show that such schemes reduce to the same IMEX R-K schemes 
for the limit equation and no parabolic stability restriction on the time step appears in the diffusive limit. 

3.1. The new approach. System p.ip can be written in the form 

u' = fl{u,v) + f2{u), 

(3.2) 

e V — g{u, V) 
where the primes denote the time derivatives and 

fi{u,v) = -{v + fip{u)x)x, h{u) = fip{u)xx, 

g{u, v) = -p{u)x -v + q{u). 

Notice that, throughout this paper, g{u, v) (and therefore gviu, u)), depends only algebraically on u, while 
it may contain differential operators acting on u. 

Now we apply an IMEX-RK scheme to system p.2[) where (/i(w, u), 0)"^ is evaluated explicitly and 
{f2{u),g{u,v))'^ implicitly. Note that if f2{u) is evaluated explicitly then by cancelation the IMEX-RK 
scheme will reduce to the typology of asymptotic preserving methods studied in [71 135) . 

In the limit e — >■ from p.2p we obtain a differential algebraic system (DAE) 

u' = fl{u,v) + /2(w), 

(3.3) 

= giu,v). 
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In order to guarantee the solvability of system p.3p we assume that the Jacobian matrix gy(u,v) is 
invertiblc, and then the DAE is said to be of index one 0. Note that if has a bounded inverse in a 
neighborhood of the exact solution, we can use the inverse function theorem in order to write 

v{t) = G{u{t)) 

for some G{u) which inserted into u' = fi{u, v) + /2(w) gives u' = G{u)) + f2{u). From now on we 
always assume that this is the case. Then, as e — > system p.2p reduces to 

u' ^h{u) + .f2{u), (3.4) 

where fi{u) = Giu)) and v = Giu). This approach will be denoted BPR approach. 

First order implicit-explicit Euler scheme that uses this approach is reported in Appendix 18.11 where 
a stability analysis is performed. In particular it is shown that as e — >■ 0, the parabolic restriction on 
time step disappears. 

In the sequel we restrict our analysis to the limit case £ — > where the main goal is to capture the 
diffusive limit. 

3.2. Analysis of TYPE A IMEX schemes. Applying an IMEX R-K scheme of type A to system 
we obtain 



Un+l = Un + At ^ Ffe) + AtY,bkf2{Uk 

k=l k=l 
s 



(3.5) 



for the numerical solution and 

fe-i 



Uk^Un + AtJ2 akj h{Uj,Vj)+ AtY^ akj /2 (C/j ) 

k 

eVfc = e^vn + AtY,akj9{Uj, Vj), 



(3.6) 



for the internal stages (notice a slight changes of notation with respect to Section ^ . 

Starting from (j3.5l) and (|3.6p . by Definition (j2.1[) and A invertible, we obtain from the second equation 
in (IXBl) 

k 

Atg{Uk,Vk) =£^Y1 ~ 
Inserting this into the numerical solution Vn+\ we make Vn+\ independent of and setting e = 0, we get 

s s 
+ AtY^ hfliUk) + At ^ bkf2{Uk) 

(3.7) 

S 

Vn+l = 'R-{(X))Vn + Aty^hkUJkjVj, 
k=l 



^The index of a DAE is the number of times one has to differentiate the funetion g to obtain a system of ODE's. For 
example, differentiating the function g, one obtains gu(u,v)u' + g^{u,v)v' = 0. If gv is invertible, system 1 13.311 ean be 
written as u' = f{u, v), v' = —g'^guf- 
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with /i(C/fe) = fi{Uk,G{Uk)), and stage values 



Uk = M„ + At ^ akj A (C/j ) + At ^ flfej /2 ([/j ) 

(3.8) 

0^g{Uk,Vk). 

ioT k — 1, . . . , s. The latter equality implies 14- = G{Uk), k ~ 1, . . . , s. 

Note that if we require that the implicit part of the scheme is stiffly accurate, i.e. if 

where = (0, . . . , 0, 1)"^, then by (EH) 

n{oo) = 1 - b^A-^l = 1 - ef 1 = 1 - 1 = 0. 

This implies that if the implicit scheme is A-stable and stiffly accurate it is also L-stable and Wn+i ~ 
V, = G([/,). 

It is interesting to note that, if we consider system (|3.ip with q{u) = 0, when e = we get a purely 
diffusive system which means that the term fi{u,v) in (|3.2p disappears. Therefore, by BPR approach, 
the IMEX R-K scheme of type A in the limit e — )• becomes a stiffly accurate DIRK scheme and hence 
no stability restriction on the time step is required in the diffusive limit, i.e. we got an unconditionally 
stable method. Another advantage of this new approach is the following. Usually the numerical solution 
(un+i,Vn+i) in (|3.7p in the case £ = will not lie on the manifold g{u,v) = since g{un+i,Vn+i) is 
not necessarily zero. But this approach guarantees that in the limit e — > we obtain a stiffly accurate 
implicit scheme and hence = Us, implying g{un+i,Vn+i) = 0. 

In the general case of systems for which q{u) ^ 0, it is fi{u, v) ^ and, by using the BPR approach, 
in the limit case £ — >■ we obtain an IMEX R-K scheme with a non vanishing explicit term in which 
the diffusion term /2(m) is treated implicitly and a classical CFL hyperbolic condition for the time step 
is required. In general g{un+i,Vn+i) 7^ even if all stage values lie on the manifold, (see the second 
equation in p.8p ). However, if the explicit scheme has the property that = Us, and the implicit 

scheme is stiffly accurate, then, in the limit as e — > the numerical solutions are projected on the manifold 
g{un+i,Vn+i) = 0, because g(u„+i,u„+i) = g{Us,Vs) = 0. 

From the above discussion it is clear that the property = Us is crucial if we want that the 

numerical solution is projected to the limit manifold as e — > 0. We emphasize that there is a class of 
s-stage explicit R-K methods for which Un+i = Us', such methods are called First Same As Last (FSAL), 
and their coefficients satisfy as,i = bi for i = 1, . . . , s — 1 and bg = 0. They have the advantage of requiring 
s — 1 function evaluations for each step (see [22] for details). FSAL methods are often used in the contest 
of embedded methods, such as the popular Dormond-Prince method (DOPRI) [18], on which MATLAB 
routine ode45 is based on. 

From the arguments above, in order to capture the limit as e ^ 0, it is important that the implicit 
part on an IMEX R-K is stiffly accurate and the explicit part is FSAL. This motivates the following 

Definition 3.1. We say that a IMEX R-K scheme is globally stiffiy accurate if b'^ = A and 
r — e'^A, with = (0, . . . , 0, 1)^, and Cg ~ Cs ~ 1, i.e. the numerical solution is identical to the last 
internal stage value of the scheme. 

From (|3.5p and (|3.6p we observe that if an IMEX R-K is globally stiffly accurate, then Un+i = Us, 
Vn+i = Vs, and therefore lim^^o g{un+i,Vn+i) = 0. 

General remarks for type A. 
• It is worth mentioning some important aspects about type A schemes. First of all, in [5] Boscarino 
emphasized that an important ingredient for the IMEX R-K schemes of type A is 6i = bi for 
all i. Such a choice provides a significant benefit for the differential component u, i.e., an order 
reduction does not appear for this component. On the other hand, conditions 

ejA^b^, ejA^b^ 

imply Oss = &s 7^ bs — Oss = which means that for a stiffly accurate IMEX R-K scheme it is 
b ^b, and therefore we expect to observe order reduction for the differential variable. 
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• It is impossible to construct a second order stiffly accurate IMEX R-K scheme of type A with 
s = 3 internal stages. The proof is given in Appendix 18.21 In practice, in order to satisfy all 
these order conditions we have to increase the number of the internal stages. In view of such 
difficulties, for type A schemes, we shall consider second order IMEX R-K schemes with s = 3 
and 6 = 6 in order to avoid the order reduction, giving up to the FSAL property of the explicit 
scheme (and with that the global stiff accuracy of the IMEX scheme) . An example is the scheme 
SSP(3,3,2) in Appendix 18.31 In this case, if q{u) ^ 0, it is g{un+i,Vn+i) 7^ as e — >■ 0. 

• Formulation p.2p in the limit case £ — yields the index- 1 DAE. Then using the same tech- 
nique adopted in [6], we can derive additional order conditions, called algebraic conditions, that 
guarantee the correct behavior of the numerical solution in the limit e — > and maintain the ac- 
curacy in time of the scheme. If the implicit scheme is stiffly accurate, such conditions becomes, 
to various order of accuracy, 

Cs = 1, (second order) , , 

e^i5== 1/2, (third order) ^ ' 

where Cg — (0, 0, 1)"^. If the IMEX schemes is globally stiffiy accurate, then p.9p are auto- 
matically satisfied, since e'^A = b^ . 

• Finally we observe that, in order to construct an order p > 3 IMEX R-K of type A and to 
maintain accuracy we have to increase the number of the classical order conditions too. Usually 
several simplifying assumptions (see [6], [8], |21| for details) could help to reduce the number of 
such conditions, but, higher orders type A schemes are more complicated to construct than CK 
or ARS schemes because of additional order conditions (see [8] ) due to the fact that c ^ c. 

3.3. Analysis of TYPE CK schemes. Similar considerations about BPR approach, explained for 
the IMEX R-K scheme of type A in the limit case e — > 0, can be reproposed here for the type CK when 
applied to the system p.Sp . with slightly modifications. Of course, if we consider the general system p.ip 
we obtain again an IMEX R-K scheme of type CK in the diffusion limit, i.e. e — >■ 0, where the diffusion 
term /2('u) is treated implicitly and a CFL hyperbolic condition for the time step is required. 

Indeed, we consider an IMEX R-K schemes of type CK where, by Definition 12. 2[ we assume that the 
submatrix A is invertible and an = 0. The Butcher tableaux of a CK scheme takes the form 












c 


a 


A 




bi 


b^ 



with a = (a2i, . . . , o.s^i)'^ and b"^ = (&2j ■ ■ ■ ,bs)- In order to simplify the analysis we consider that the 
implicit part of the scheme is stiffly accurate. Under this circumstance it is easy to prove that 



bi+iFa = 0, (3.10) 

where a = —A~^a (see [5] for details). 

Then, considering a scheme of the type CK, the second equation in p.6p becomes 



k 

e^Vk = e^v,, + Atakig{un,v,,) + AtY,akj9{Uj,Vj). (3.11) 

i=2 

with k = 2, s. 

Now multiplying by LOkj, where LOkj sue the elements of the inverse of A, and summing on /c, we 
obtain 

Atg{Uk,Vk) ^ e'^^UkjiVj ~ Vn) + Atakg{un,Vn), for fc = 2,...,.s 

J=2 



where 



1=2 1=2 

Insertmg the expression Atg{Uk, Vk) into the second equation in p.Sp we obtain 

e^u„+i = e^7^(oo)^;„ + At^bkUJk]V.j + Ai 61 + ^ 6^0^ j g{u„, v„) (3.12) 

A:=2 \ k=2 ) 

Then by p.lOp the last term in the second equation in p.l2p drops and in the hmit case for e = we 
can write 

s 

Vn+\ = TZ{(X))Vn + Aty^^bkUJkjVj, 

k=2 

with 

giUk,Vk) = akgiun,Vn), for fc = 2,...,s. (3.13) 

Note that, for IMEX R-K schemes of type CK, the stabihty function TZ{z) of the imphcit part of the 
scheme takes the form 

TZ{z) = I + z{bi + - zA)-\ls-i + za)) 

(3.14) 

= (61 - b^A-^a)z + (1 - b^A-^1,^1 + b^A-^a) + 0{-). 

z 

We obtained this resuh, by applying one step of the imphcit part of the scheme to the test problem 
y' = Ay, y(to) = 1, with A G C and 1^ = (1, . . . , 1)^ S W. 

Thus, the only stiffly accurate condition, i.e. eJ_iA = b^ is not enough to guarantee that limz^ooTi-iz) = 
and then an additional condition is required for the implicit part of the scheme, (for details sec [5]). 
This is expressed by the following 

Proposition 3.2. // 

— e^-iA~^a — "^^uisjaji — 0, (3.15) 
then 7^(oo) = 0, where e^-i = (0, . . . , 0, 1)"^ € . 

Proof. In fact, assuming A invertible, we get h^A^^ = ef_]^ and when z — 00, from p.l4p we obtain 
7?.(oo) = b^A~^a = —e^_iA~^a, which is zeros if p.isp is satisfied. □ 

Note that the previous Lemma implies that ag ~ —e'^_iA~^a — and by p.l3p with k — s wc obtain 
g{Us, Vs) = 0, then the last stage values lie on the manifold g{u, u) = as e 0. Now wc observe that if 
the IMEX R-K scheme of type CK is globally stiffly accurate, we obtain from (|3.12p and (|3.7p w„+i = Us 
and Vn+i = Vs and therefore (7(w„+i, w„+i) = with w„+i = G'(u„+i). 

Since an IMEX R-K schemes of type ARS is a particular case of the type CK where the vector a = 0, 
then the same results hold true. 
General remarks for type CK. 
• IMEX CK schemes [11] are attractive because of their good properties. The implicit part of this 
scheme is singly diagonally implicit Runge-Kutta (SDIRK) with Oii = 7 > for i = 2, ...,s and 
differs from the classical SDIRK one because an = 0. In [11] such implicit schemes are called 
explicit singly diagonally implicit (ESDIRK). A consequence to set an = is the possibility to 
guarantee stage-order q higher than the in the case of SDIRK, for which q — \. Moreover here wc 
consider schemes that arc stiffly accurate according to Definition 13.11 Such schemes will project 
the solution on the manyfold in the limit of infinite stiffness. For these schemes b ^ b, so one 
of the so-called simplifying conditions cannot be applied [8]. Here we require that q — Ci for all 
i = 2, s; this choice will reduce the number of coupled order conditions. 
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4. IMEX-Finite Difference schemes. When constructing numerical schemes, one has also to 
take a great care in order to avoid spurious numerical oscillations arising near discontinuities of the 
solution. This is avoided by a suitable choice of space discretization. To this aim it is necessary to use 
non-oscillatory interpolating algorithms, in order to prevent the onset of spurious oscillations (like ENO 
and WENO methods), see |40]. Moreover the choice of the space discretization may be relevant for a 
correct treatment of the boundary conditions. 

In this section we emphasize some requirements about the space discretization of the system 
We remark that the dissipative nature of upwind schemes [Ml |3S] depends essentially on the fact that 
the characteristic speeds of the hyperbolic part are proportional to 1/e. On the other hand central 
differences schemes avoid excessive dissipation but when e is not small or when the limiting equations 
contain advection terms may lead to unstable discretizations. In order to overcome these well-known facts 
and to have the correct asymptotic behavior we fix some general requirements for the space discretization. 

1. Correct diffusion limit. Let us consider system p.ip with q{u) = 0. In the limit case e — the 
therm v + dxp{u) from the second equation. If we want that v + fi{e)dxp{u) — >■ also in the 
first equation, we need to use the same space discretization for the term dxp{u) and require that 
/.(O) = 1. 

2. Compact stencil. Among the advantages of our approach there is the possibility to have a scheme 
with a compact stencil in the diffusion limit e — > 0. This property is satisfied if point 1) is satisfied 
and we use a suitable discretization for the second order derivative that characterize the diffusion 
limit. 

3. Shock capturing. The schemes when q{u) ^ should be based on shock capturing high order 
fluxes for the convection part. This is necessary not only for large values of e but also when 
we consider convection-diffusion type limit equations with small diffusion. The high order fluxes 
are then necessary for all space derivatives except for the second order term fi{e)dxxP{u) on the 
right-hand side. 

4. Avoid solving nonlinear algebraic equations. In order to achieve this the implicit space derivative 
dxP{u) in the second equation must be evaluated using only nodal values of u which can be 
obtained from the solution of the flrst equation. 

The above properties are satisfied for example using high accuracy in space obtained by finite difference 
discretization with Weighted-Essentially Non Oscillatory (WENO) reconstruction, [40] . 
System (j3.ip may be written in the form 

Ut + {v + pp{u)x)x = pl-P{u)xx, 

(4.1) 

^* = ^ - + P{u)^)) ■ 

with ji = iJL{e) introduced in Section 2. The terms on the right-hand side will be treated implicitly. For 
large value of e the explicit flux is just (w,0)'^, while for small values of e it is {v + p{u)xtQ)'^ ■ Here we 
describe a finite difference WENO scheme for a system of the form 

Ut + F{U)x = R{U), 

and apply it to the system (|4.1[) with 

F{U) ^ {v + ^Jip{u)x,Qf , 
R{U) = {pLp{u)xx, \ {q{u) - {v+p{u)x))). 
As e — > oo and — > 0, the system becomes 

Ut+Vx= 0, 

and the characteristic speed of the system is A ~ (twice). As e — !• and 1, w + ^p{u)x q{u) and 
the system relaxes to the equation 



ut + q{u)x = p{u)xx 
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and the characteristic speed of the left hand side is given by A = q'{u). 

Conservative finite difference for system (|4.ip are written as follows, [35] 



where Uj{t) « U{xj, t) is an approximation of the pointwise value of U at grid nodes, and the numerical 
flux at cell edge x^j^i is computed as follows 

The function F^{x) and F~^-^{x) are suitable reconstructions defined, respectively, in cell j and in cell 
j + 1. They are obtained as follows. First, we assume that the flux can be split into a positive and 
negative component 

F{U)^ F+{U) + F-{U), 

with X(yuF+{U)) > 0, A(V[/F-(C/)) < 0. The quantity F^ F'^'^Uj) are computed at ceU center. 

Then Ff{x) are reconstructed from {F^} using high order essentially non oscillatory reconstruction, 
such as ENO or WENO, that allows pointwise reconstruction of a function from its cell averages, (see, 
e.g. [40j for details). 

The flux F may contain derivatives. For example the flrst equation in system ()4.ip contains p{u)x. 
Such terms are computed by point-wise WENO reconstruction. 

In all our examples wc used the simple local Lax- Friedrix flux decomposition, i.e. F^{U) ~ ^{F{U) + 
aU), F'{U) ^ \{F{U) - aU), a > maxy |p(V[/F)|, VA e R"^™, where p{A) = maxi<,<„ |A,(A)| 
denotes the spectrum radius of matrix A, and the max defining a is taken for U varying in a suitable 
range in a neighborhood of each cell. In our test case we chose a — 1, since as e — >■ cx), p{W{jF) — and 
in our numerical test q{u) is either 0, u, or with U ranging in [0, 1], therefore < 1- 

We remark here that the choice of a is based on the characteristic speeds of the limit convection- 
diffusion equation, while a more detailed analysis is needed to justify its use in intermediate regions, 
for which the characteristic speeds can be much higher, and the stabilization that compensates for the 
apparent violation of the hyperbolic CFL condition comes from the implicit treatment of the diffusion 
term. 

Furthermore for large value of e, (e.g., e = 1), we want to avoid adding and subtracting terms which 
may cause loss of accuracy. For a semidiscrete scheme the function fi will depend also on the grid space 
Ax. A simple choice for /x is given by 

fi = cxp(— e^/Acc) 

which is what we used in all our numerical tests. 

For the diffusion term p{u)^x we used the standard 2-nd order finite difference technique for second 
order time discretization, and the standard 4-th order finite difference technique where 3-rd order time 
discretization are used. 

5. Numerical examples. In this section we test several second and third order IMEX R-K schemes 
presented in the literature that satisfy the algebraic order conditions (|3.9p and conditions in Definition 
3. Usually, IMEX time discretization are identified by an acronym (e.g. the initials of the authors), and 
three numbers {cfe, ct/, p) denoting, respectively, the effective number of stages (in practice the number 
of function evaluations) of the explicit and implicit scheme and the classical order of accuracy. 

Below we list the IMEX R-K schemes used in the numerical tests. 

• SSP(3,3,2): derived by Pareschi, Russo [Ml, it is a second order IMER R-K of type A, the explicit 
part is strongly stability preserving, while the implicit part is stiffly accurate. In accordance with 
the proposition 18.11 in the Appendix 18.21 this scheme is not globally stiffly accurate according to 
Definition 13.11 

• ARS(2,2,2): derived by Asher, Ruuth, Spiteri [1], it is a second order scheme, the double Butcher 
tableau of this scheme is reproduced in Appendix 18.31 Note that this scheme is globally stiffly 
accurate according to the Definition 13.11 and satisfies the additional order conditions p.9p . 
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• ARS(4,4,3): derived by Asher, Ruuth, Spiteri [T], it is a third order scheme, the double Butcher 
tableau of this scheme is reproduced in Appendix 18.31 Similarly to ARS(2,2,2) this scheme is 
globally stiffly accurate according to Definition 13. f I and satisfies the additional order conditions 

• BPR(3,5,3) introduced in this paper is a third order IMEX R-K scheme of type CK and globally 
stifHy accurate according to Definition l3.f I This scheme has s = 5 internal stages, ctb = 3 explicit 
stages and aj = h implicit stages. The additional order conditions p.9|) are satisfied. This scheme 
is more efficient than ARS(4,4,3) for the explicit part, but less efficient for the implicit one. In 
many cases the computation of the explicit term is more expensive than the solution of the 
implicit step, thus resulting in an overall improvement in efficiency per time step. Furthermore, 
the coefficients of the explicit scheme are all non negative, which is an advantage for the stability 
properties of the scheme. We reproduced the coefficients of this scheme in Appendix 18.31 

In all the computations presented in this paper we denote each scheme with an acronym indicating the 
IMEX scheme and the type of space discretization. 

Space discretization is identified by a short name containing the order of accuracy in space; for 
example, WEN053 (or WEN032), see for details [40], means a fifth (or a third) order reconstruction 
which reduced to third (or second) order near singularities and CdS2 stands for second order central 
discretization scheme. 

We remark here that all the analysis performed in the paper and the numerical tests are performed 
under the assumption that the initial data is well-prepared, which means that the initial condition lies in 
the limit manyfold as £ — > 0. If this condition is not satisfied, then a loss of accuracy is observed, unless 
some initial layer fix is adopted. Schemes of type A are more robust against this problem, as is described 
in [36] . 

5.1. Convergence test. In this section we investigate numerically the convergence rate of the 
second and third IMEX R-K schemes introduced before for a wide range of the parameter e. To this aim 
we apply these schemes to simple prototype hyperbolic system p.ip . with initial conditions chosen in 
such a way that the exact solutions is smooth and docs not present a rapidly varying transient for small 
values of e. This is achieved in practice by imposing that the initial condition satisfies the limit relation 
between u and w as £ — )- 0. 

Numerical convergence rate is calculated by the formula 

V = log2(-BAti/-BAt2), 

where -Ba*! and £'At2 are the global errors computed with step Ati = 0(Ax), and At2 = Ati/2. In the 
following tests we put = 10^^ and we choose Ai w Ax. 

For the first test we set p{u) — u and q{u) = 0. Then we get 

— ^'^xx ~t- fiUxx: ('5 1') 

e'^vt = ~Ux -V, 

that in the limit ease, £ = and /i = 1 leads to the linear diffusive problem 

ut=Uxx, u{x,0) ^ uo{x). (5.2) 

We use periodic boundary conditions with uq{x) ~ cos(a;), and x £ [0 , 27r], so that u{x, t) ~ uq{x) exp(— 
is an exact solution of (|5.2p . The final time is T = I and At = 0.5 Ax. 

The results are reported in Table 5.1 and 5.2 showing that the expected convergence rates are reached 
for the w-component. 

Next we set p{u) — q{u) — u and consider the following system 

Ut+Vx= fiUxx - fiUxx 

e^vt + Ux = -{v - u), 

where the limiting behavior is given by an advection-diffusion equation. We use periodic boundary 
conditions with the initial data u(x, 0) = exp(— (I + cos(x — Tr))/a), v(x, 0) = u{x, 0)(1 — /zsin(a; — 7r)/cr) 
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Table 5.1 
Convergence rate for 11 tn Licio-TiOTTn, 



ARS(2,2,2)-CdS2 
N Laa{u) Order 


SSP2(3,3,2)-CdS2 
N Lac(u) Order 


ARS(2,2,2)-WEN032 
e^A = IT 
N Loo(u) Order 


SSP2(3.3,2)-WEN032 
A ^ 6^ 
N Lao(u) Order 


20 7.800e-03 
40 1.873C-04 2.05 
80 4.597e-04 2.02 
160 1.138e-04 2.01 
320 2.833e-05 2.00 


20 2.906C-02 
40 7.979e-03 1.86 
80 2.039e-03 1.96 
160 5.120e-04 1.99 
320 1.274e-04 2.00 


20 4.820e-03 
40 1.492e-03 1.69 
80 4.124e-04 1.85 
160 1.082e-04 1.93 
320 2.760e-05 1.97 


20 4.697e-03 
40 1.483e-03 1.66 
80 4.102e-04 1.85 
160 1.074e-04 1.93 
320 2.748e-05 1.96 



Table 5.2 
Convergence rate for 



ARS(4,4,3)-WEN053 
ejA = IT 
N Lac(u) Order 


BPR(5,5,3)-WEN053 
A = bT 
N Loo{u) Order 


20 1.810e-02 
40 3.365e-03 2.42 
80 5.349e-04 2.65 
160 5.960e-05 3.16 
320 5.968e-06 3.31 


20 1.639C-02 
40 3.099C-03 2.40 
80 5.167e-04 2.58 
160 5.821e-05 3.14 
320 5.949e-06 3.29 



with <T = 0.05 and /i = 1, on the spatial interval [0, 2tt], at the final time T = 0.3 anfd At = O.SAx. As 
reference solution we use the truncated Fourier representation of the exact solution 

k— — oo k—--oc 

with Uk{t) and Vk{t) satisfying 

ilk = -ifcVfc, , . 

e^Vk^~ikUk + Uk-Vk. ^ ■ ' 

For each k, system (|5.3p can be written as a 2 x 2 constant coefficient homogeneous system which 
can be solved exactly. The results are given in Table 5.3 showing that again the expected convergence 
rates are reached for the u-componcnt by all schemes. 

The above convergence analysis has been performed in the limit e — > 0, therefore we might expect 
a degradation of the accuracy for intermediate regimes as in the case of hyperbolic relaxation when 
the classical order is greater then two [SJ [51 [TT]. Furthermore, from the practical point of view, the 
understanding of this phenomenon is essential in situations where one is interested in the construction of 
higher order methods. 

Figure 15.11 shows the convergence rates as a function of using different values of ranging from 
10"*^ to 1 and At ~ Ax. Second order schemes ARS(2,2,2)-CdS2, SSP2(3,3,2)-CdS2 have the prescribed 
order of accuracy uniformly in (upper left panel). Instead, ARS(2,2,2)-WEN032 and SSP2(3,3,2)- 
WEN032 present a degradation of accuracy at intermediate regimes (upper right panel) . 

A similar lack of convergence in intermediate regimes is observed for both the third order schemes 
ARS(4,4,3)-WEN053 and BPR(3,5,3)-WEN053 (lower left panel). This resuhs have a very different 
nature than the accuracy degradation observed in IMEX schemes applied to hyperbolic systems with 
stiff relaxation [5]. A plausible reason here appears to be that in intermediate regimes the two terms 
which are added and subtracted in the equations, i.e. ±fj,p(u)xx, arc discretized in two very different 
ways: one is computed inside the flux (ARS(4,4,3)-WEN053 and BPR(3,5,3)-WEN053), and the other 
one (ARS(4,4,3)-WEN053* and BPR(3,5,3)-WEN053*) is discretized by a discrete one dimensional 
Laplacian, therefore the two terms do not almost cancel each other. Although in certain regimes such 
problem could be solved by treating the term p{u)x out of the flux (see, for example, the result in the 
lower right panel of Figure [5?T] for BPR(3,5,3)) and discretizing both terms ±^p{u)xx in the same way, 
this may compromise the cancelation of the quantity q{u) — v — p{u)x in the flux. A general understanding 
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Table 5.3 

Convergence rate for u in Laa-norm in the convection- diffusion limit. 



ARS(2,2,2)-CdS2 
N Loo(u) Order 


SSP2(3.3,2)-CdS2 
N Loo(u) Order 


ARS(4,4,3)-WEN053 
ejA = bT 
N Loa(u) Order 


BPR(5,5,3)-WEN053 
N Loo(u) Order 


40 3.867C-03 
80 9.457e-04 2.03 
160 2.330e04 2.02 
320 5.798e-05 2.00 


40 2.615e-03 

80 6.243e-04 2.06 
160 1.543e-04 2.01 
320 3.850e-05 2.00 


40 4.297e-04 

80 5.770e-05 2.89 
160 7.922e-06 2.86 
320 1.256e-06 2.65 


40 8.300e-04 

80 1.167e-04 2.83 
160 1.603e-05 2.86 
320 2.230e-06 2.85 



and a robust treatment of intermediate regimes is beyond the scope of the present paper and requires 
further investigation. 




Fig. 5.1. Convergence rate in Lao-norm versus for various schemes. Lack of convergence for intermediate values 
of e is evident in the upper right and louier left panels. The lovjer right panel shows results obtained by a scheme in ujhich 
the explicit and implicit term iip(u)xx is discretized spatially in identical way, namely by a second order discrete Laplacian. 

5.2. Shock test cases. In this section we apply the scheme to problems with discontinuous initial 
data, that in the limit as e — > reduce to convection-diffusion equation. Notice that in the relaxed limit 
the scheme becomes an IMEX scheme for the limit equation. This test is used to check both the shock 
capturing properties of the scheme, and its relaxation to an IMEX scheme for the limit equation. 

In the rest of the section we will consider the third order BPR(3,5,3)-WEN053 scheme. 
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Fig. 5.2. Numerical solutions at time t = 0.25 in the rarefied regime (e^ = 0.7) with At = 0.5 Ax and Ax = 0.01. On 
the left-hand side the mass density u (a) and on the right-hand side the flow v (b). Solid line is the reference solution. 



First we consider a purely diffusive linear problem. We solve a Riemann problem, in the rarefied and 
diffusive regime for system 

e^vt + Ux = -V. 

We take the following initial data 

UL = 2.0 VL =0, -Kx < 0, 
ur = 1.0 vr = 0, < .T < 1. 

As e goes to zero we get ut = Uxx, i-e. the problem becomes a classical Riemann problem for the heat 
equation. 

In order to test our scheme we compute the numerical solution in the rarefied (e^ > Ax) regime and 
in the diffusive (e^ ^ Ax) regime. This means that when is very large (i.e., rarefied regime) /i is very 
small, and on the other hand when e is very small (i.e., diffusive regime) fi is equal to 1. 

We set = 0.7 for the rarefied regime and — 10~® for the diffusive regime (or stiff regime). The 
numerical solution for u and v in the rarefied (Fig l5.2|) and diffusive regime (Fig l5.3|) are depicted with a 
reference solution obtained using a fine spatial grid of iV 2000 cells. 

As boundary conditions we set w = Oata; = ±l. Compatibility with the system gives Ux — at 
a; = ±1. Notice that the characteristic variables for this problem are S,± = u ± sv, therefore condition 
w = Oata:; = ±lis equivalent to impose ^+ = ^_ at the boundary. For such a reason we denote these 
boundary conditions as reflecting boundary conditions. 

The solution is reported at final time t — 0.25 in the rarefied regime (Fig. 15. 2p and t — 0.04 in the 
diffusive regime (Fig. 15. 3|) . In the figures we observe that the scheme captures well the correct behavior 
of the solutions both in rarefied regime where it provides an accurate description of the shock without 
oscillations near the discontinuities, and in the diffusive regime where the numerical solution matches 
accurately the reference solution. 

Finally we consider the nonlinear Ruijgrok-Wu model, [38], (for details see |25j ) 

/ \ '^XX / \ '^xx 

2ko 2ko ^^^^^ 



e vt+ Ux = -2ko 



C , n 2 2 \ 

V — —(u ~ e V ) 
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Fig. 5.4. Numerical solutions at time t = 2.0, (a) in the intermediate regime (e^ = OA) and (b) parabolic regime 
(e^ = 10~^) with At = 0.25Aa; and Ax = 0.2. Solid line is the exact solution. 



where we added and subtracted to the original model the quantity fiu^x on the right-hand side of the 
first equation and C is a constant, (in our test we chose C — 1). In the diffusive limit e ^ 0, the second 
equation provides 



1 2 1 

V — -U — —T—Ux, 

2 2kn 



and we get the limiting viscous Burgers equation 



Note that through BPR approach the IMEX R-K scheme in the diffusive limit relaxes to the same 
IMEX R-K one for equation ()5.5p where the convection term is treated explicitly and the diffusion term 
implicitly. 

The exact solution to the shock- wave problem has been given in |38| . The initial conditions are two 
local Maxwellians characterized by 

UL=2.0, -10<a:<0, 
ur = I.O, < X < 10, 

with w = [(I + M2e2)i/2 _ l] /£2^ 

In Figure WM we show the computed solution for the mass density u in the intermediate (e^ = 0.4) 
and parabolic (e^ — 10~^) regimes versus the exact solution. As it can be seen once again, the scheme 
gives an accurate description of the viscous shock profiles. 
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6. Application to transport equations. In this section we apply the IMEX schemes derived in 
the first part of the manuscript to the case of neutron transport equations [55J [5^1 121] • 

We consider the multidimensional transport equation under the diffusive scaling. Let /(t,x, v) be 
the probability density distribution for particles at space point x € M'', time t, traveling with velocity 
V S fi C M'' with d-v ~ S. Here il is symmetric in v, meaning that g{v)dv = for any function g 
odd in V. Then / solves the non dimensional linear transport equation 



(6.1) 



where a = ^(a;) is the total cross section, as ~ crs{x) is the scattering coefficient. Here Q = Q{x) is a 
source term and e the mean free path. Typically, cr^ = a — e^OA where a a = <yA{x) is the absorption 
coefficient. Such an equation arises in neutron transport [12] . radiative transfer jl4j and wave propagation 
in random media |39| . etc. In all these applications, the scaling appearing in (j6.ip is typical, and gives 
rise to a diffusion equation as e of the form [55] 

ftp = ^ y" V • V.^, ■ Vxp^ dw - GAP + Q, 

where p = (l/S) J^^ f dv. We refer, for example, to [HIJ] and the references therein for rigorous mathe- 
matical results concerning diffusion limits of transport equations. 

6.1. Problem reformulation. Consider now the one-dimensional transport equation 

edtf + vd^f = 1{yJ, " '^•^) + ^^-^^ 



with XL < X < xr and boundary conditions 



f{t,XL,v)^FL{v), for w > 0, 
f{t,xn,-v) ^ Fr{v), for v > 0. 



(6.3) 



In |26| the authors proposed a method based on the even-odd decomposition f = r + ej where r = 
i(/(u) -I- /(—«)) and j = '^{fW) — fi—v)), that splits the equation (|6.2p as two equations, each for v > 

1 /cr« 



edtfiv) + vd^fiv) ^-l^-J- J ^ fdv' ~ <jf{v)j + eQ, (6.4) 
edtf{-v) - vd,f{-v) = J ( Y fdv' - af{-v)^ + eQ. 



Adding and subtracting these two equations leads to 



dtr + vdxj^ ^{r - p) - aAT + Q, (6.5) 



where 



As e — 0, system (|6.5p gives 



dtj + -^OxT = 5-J - (TAJ, 



fi 

P = I rdv. (6.6) 



r = P, j = -{v/a)dxr. 

Applying this to the first equation of (|6.5p and integrating over v we get the diffusion equation 



9tP = ^dxxP - (7AP + Q- (6.7) 
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To get boundary conditions for r and j we use relations 



r + £.?U=xi = Fl{v), r ~ ej\^=.^„ = Fr{v). (6.8) 

For arbitrary value of e, the compatibility conditions imposed at the boundary is quite complicated. 
However, for small values of e, which is the case we are interested in, the treatment simplifies, because 
when £ — >■ 0, j = —{v/a)dxr, then applying this in (|6.8p one gets 



r 



svd^rl^^x^ = Fl{v), r + evd^rlx^^^ = Fji{v). (6.9) 



Such boundary conditions will avoid a boundary layer in the limit case e — > 0, therefore the numerical 
boundary conditions are obtained by discretizing Eq. (|6.9p . As done in [55], the boundary conditions 
have been applied using a second order implementation of equation (|6.8p based on central differences. 
Extensions to higher-order implementation of equations (|6.9p and to different boundary conditions are 
not considered here and will be investigated in a forthcoming work. According to this, for our tests we 
chose a second order scheme in space and time, because the boundary conditions are discretized to second 
order accuracy. 

Now we start from system (|6.5|) and, adding and subtracting the quantity v^dxxf/cr in the first 
equation, we reformulate the problem in the equivalent form 



Otr = -vdx{j H ) 5-(r - p) - <JAr + Q + ^i{e)v^ 



Explicit Implicit to -i r\\ 



Implicit 

The new system is then discretized in time using an IMEX R-K scheme as described in the first part of 
the manuscript. The implicit-explicit integration process is emphasized in (|6.10p . 

We remark that this new formulation of the diffusive relaxation system (|6.5p is such that when e 
tends to zero the system relaxes towards (|6.7p . From a numerical point of view the new formulation has 
several advantages. In particular, as e — >■ 0, the IMEX R-K scheme applied to system (|6.10p originates a 
fully implicit scheme for solving the diffusion equation (j6.7p . 

In our numerical tests in order to obtain uniformly accurate second order scheme both in space and in 
time for the BPR approach, we consider SSP2(3,3,2)-WEN032. We remark that we can obtain analogous 
results considering ARS(2,2,2)-WEN032. 

The equations are discretized in space and velocity, i.e. r(xi,Vmitn) ~ ffm where {vm} are chosen 
to be the N.^, positive nodes of the Gauss-Legendre quadrature formula, with 2Ny nodes in the interval 
[— 1, 1] while Xi = Ax{i — l/2) for i = 1, Np. Note that the computation of the k-th stage of the implicit 

equation r^^^, requires the quantity p^'^'' in the implicit part in (|6.10p . Such quantities are obtained as 

follows. Assume we have computed r^^'^ for / — 1, fc — 1, then r^'^i is obtained from 



(k) ^ 



t:''+^i-^^i^ir^.-fr^)--Ar^. + Q) (6-11) 



discretizing system (j6.10p without the quantity ^{e)v'^ Oxx^ / <J in the implicit and explicit part. The 
quantity ^ represents the contribution of the first fc — 1 stages. Then, in order to compute p] we 
apply Gauss quadrature on both sides of ()6.1ip (i.e. multiply by the weights Wm and sum over m), setting 

(k) 

= 0. In this way we obtain an equation for p]^ that can be explicitly solved, and such value is plugged 
in (|6.10l) in order to compute r^^^. 

6.2. Numerical results. In this section we shall consider some transport problems in slab geometry. 
We will present the transient and the steady state solutions. We remark that in all the test problems we 
have used Ny = 8 thus the standard 16 points Gaussian quadrature set for the velocity space. In all the 
tests the initial distribution is /(x,v,t = 0) = 0. 

We emphasize that, besides uniform accuracy in e, our approach allows to choose larger time steps, 
since there is no stability restriction on the time step. As we will show, this permits to obtain numerical 
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results at a lower computational cost compared to other approaches presented in the literature that lead 
to explicit schemes for the underlying diffusion limit with a parabolic CFL stability restriction [55] . 
Nevertheless, in order to get an accurate resolution of the behavior of the solution, smaller time step may 
be necessary. 

Depending on the regime of the parameter e, we compare numerical solutions to a direct implicit 
discretization of the diffusion limit (|6.7|) when e tends to zero, whereas for intermediate values of the 
parameter e we compute a reference solution using a much finer grid in space. 

In the next tests we compare the results obtained by the new approach versus the results given by 
Jin et al., in [2^, here denoted by JPT. We refer to [26l [28l |32] for similar results where the limiting 
scheme is explicit and so in diffusive regions requires At (Ax)^. In all figures we use notations A^^ and 
Np to denote the number of time steps and grid points in space respectively. 

Problem I:. 

xe[0,l], Fi(v) = l, Fr{v)^0, 

as = l, <7A^0, Q = 0, £ = 10-8. 

The numerical results are reported in Figures |6TT] and [6?2] fa) at different times t = 0.01, 0.05, 0.15 with 
Np = 40, and at t = 2 with Np = 20 where the steady state is reached. In this problem we see that in 
both cases, the results in the transient and steady state solutions show a good behavior with the correct 
diffusion limit. The exact diffusive solution has been computed by (|6.7p with Np ~ 200. As expected 
both JPT and BPR results are very close to the exact diffusive solution at any times. Note however 
that thanks to the better stability properties in this regime BPR scheme is about 4 times faster then the 
explicit method. 

Problem II:. This is a two-material problem used in |26 1 I32 [ [28] where in the purely absorbing region 
[0, 1] the solution decays exponentially whereas in the purely scattering region [1, 11] the solution is 
diffusive, the parameters are the following 

a; e [0,11], Fl{v) = 5, Fr{v) = 0, 
erg = 0, CTA 1, Q = 0, e = 1, for x G [0, 1], 
CTs = 1, CTA=0, Q = 0, £ = 0.01, for [1,11]. 

An interface layer is produced between the pure absorbing region and the scattering one. Two meshes 
are used in the domain [0, 11], a thin mesh Ax = 0.05 in [0, 1] and a coarse mesh Aa; = 1 in [1, 11], 
which means that between the interface layer we have to use a space discretization with a non uniform 
mesh. Finite volume, rather than finite difference, is used in this case. Since we restrict to second 
order accuracy, the point-wise values of the source is identified with its cell average. The high order non 
oscillatory reconstruction is performed by a WENO approach for non uniform mesh. For our numerical 
tests we used WEN032 reconstruction. 

At time t = 150 the solution has reached the steady state and the results are presented in Figure |6^ 
We computed the reference solution with a very fine discretization Np — 400 using a uniform mesh in all 
the domain [0, 11]. The numerical schemes provide a good description for the solution in the absorption 
and diffusive regions, in fact, we observe that for the steady state. Figure [Ql (b) . JPT and BPR results 
are close o the reference solution, except at the interface where a slight difference is observed. In this 
case BPR scheme is about twice times faster then JPT approach. 

Problem III:. Concerning this problem we present two different situations (see [26l[32]) with non- 
isotropic boundary conditions that generate a boundary layer: 

a;e[0,l], Fl{v)=v, Fr{v)=^0, 
as = 1, (TA = 0, Q = 0, £ = 10"^ 

First in an intermediate regime with e = 10^^ and then in a more diffusive regime with e = 10-**. Using a 
coarse discretization Np ~ 25 the boundary layer is not resolved, but we observe that the two approaches 
accurately capture the solution inside the domain (in Figure ()6.3p we have restricted the numerical and 
the reference solution to the interval [0,0.5]). The reference solution has been obtained using a fine 
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discretization Np = 400 and the boundary layer is resolved. The results are plotted at time t = 0.4 in 
figure [^751 The higher efficiency of the present method results in an improved time step ratio of a factor 
4. 



Fig. 6.2. The steady state solution for the Problem I and II. 




(a) Problem I. The steady state solution of the mass den- 
sity p. JPT (o), At = 0.001, Ns = 2000. BPR (0), 
At = 0.0035, Ns = 561. The exact diffusive solution is 
represented by the solid line. 



(b) Problem II. The steady state solution of the mass 
density p with e = 1 and Ax = 0.05 on (0, 1) and with 
e = 0.01 and Aa; = 1 on (1,11). JPT (o), with At = 
0.025 and = 6000. BPR (0), with At = Ax, with 
Ns = 3000 



7. Conclusions. In this manuscript we have presented a general way to tackle diffusion limit for 
hyperbolic and kinetic problems which permits to obtain accurate and efficient schemes both in rarefied 
and diffusive regimes. The new approach, in particular, give rise to a fully implicit method for the 
diffusion component of the limiting system. This is obtained without solving nonlinear systems of implicit 
equations but by a suitable blending into the IMEX R-K method of a fully implicit solver for the limiting 
diffusive system. Numerical results show that this approach is able to capture the correct asymptotic 
behavior of the system at a lower computational cost compared to other approaches that lead to explicit 
schemes for the underlying diffusion limit, because we removed the parabolic CFL restriction, common 
to most approaches in the literature. 
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Fig. 6.3. Problem III. The mass density p. JPT (o), with e = 10"^, Ax = 0.04, At = 0.001, Nj, = 25, Ns = 400. 
BPR (()), Az = 0.04, At = 0.002, Np = 25, A^^ = 100. Left e = 10"^, right e = lO"' The reference solutions are 
represented by the dash dot line. 



The method here presented is based on the use of IMEX Rungc-Kutta methods, however exten- 
sion to more general additive Runge-Kutta schemes are naturally possible. In fact from our problem 
reformulation 



explicit implicit 



implicit 



(7.1) 



we can clearly combine different implicit solvers to tackle the "highly" stiff component g{u,v) which 
originates the algebraic condition g{u,v) = giving rise to the equilibrium projection v = G{u) and the 
"mildly" stiff component f2{u) corresponding to the limiting diffusive term. We leave this possibility to 
future research directions. 

8. Appendix. 

8.1. Stability analysis of first order IMEX schemes. For the subsequent analysis we restrict 
to the linear case p{u) = u and q{u) — 



Ut = -{V + ^iUx)x + ^J.U^ 



(8.1) 



We now look for a Fourier solution of the form u ~ u{t) exp(z^x), v = v(t) exjp{i^x) and inserting the 
ansatz into systems (jl.I[) , the evolution equations are 



e^vt 



(8.2) 



It is convenient to rewrite the system using the variable w ~ ^iv/S, in place of v so the system becomes 

Ut = C^(l& + /iu) — S.^fJ'U, 



e^Wt = —u — w. 

We apply the first order IMEX method based on explicit and implicit Euler schemes to get 
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(8.3) 



(8.4) 



which after manipulation can be written exphcitly in the form 



1 + At^^^l ' 



(8.5) 



At 



-w 



(e2 + Ai)(l + AtCV) e^ + At ' 

In order to study the stabihty of the method we compute the eigenvalues of the stability matrix 

/ 1 A^g^ \ 



R = 



1 

At 



V 



1 + At^^^l 

s^i+Atef^)-At^e 



(8.6) 



We obtain the expressions 
A± 

with 



'e^ + At (£2+At)(l + At$V) / 
i {l + a(l + /3) - /3 ± V[l + a(l+/?)-/3]2-4a} 



•7) 



Ate 



e2 + At' 



1 + /iAt^2 ■ 



It can be shown that |A±| < 1 when 



a < 



1-2/3 + /?^ 



1 + 2/? + /32 ■ 

The above inequality involves a third order polynomial in At 

- eit^ ~ l?At^ + 2^2(2e2^2^ ^ ^ _ ^^^^2 ^ ^4^2^2 _ i) < 0. (8.8) 

The roots of this polynomial are given by Tq = and 

Condition (|8.8p can be satisfied only if the last two roots are positive. This is guaranteed when 2e|^| < 1 
and so we have the time step restriction At < r_. The largest stability region is obtained when jj ~ I 
for which we get 



2^2\ 



CAt < 



1 (1 - A^e 
4 i2^ 



(8.9) 



This relation is interpreted as follows. For a fixed ej^j < 1/2, the restriction on the time step is of 
parabolic type, since |^| ~ 1/Ax is the maximum Fourier mode represented on a grid of spacing Ax. The 
restriction on At/ Ax^ is less and less severe as e — )- 0. 

The implicit treatment of the second equation stabilizes the explicit treatment of the first one, 
provided e is sufficiently small. 

8.2. Analysis of second order stiffly accurate schemes. We have the following result: 
Theorem 8.1. Consider an IMEX Runge-Kutta scheme of type A. Then there exist no second-order 

tree stage scheme satisfying the conditions h"^ = and = ejA. 

Proof. We consider the classical second order conditions 



b^e = 1, b^e = 1, 
6^5=1/2, 6^c=l/2, 
6^c=l/2, 6^5=1/2, 
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(8.10) 



with c = Al and c = Al and the conditions b'^A ^ = and = 

For s = 3 the Butcher tableau of a stiffly accurate IMEX R-K of type A is 















Cl 


Cl 










C2 








C2 


C2 - a22 


022 





1 




b2 





1 




&2 


7 




hi 


b2 







bi 


&2 


7 



(note that stiff accuracy imphes 
exphcitly written 



1) and from (jS.lOp the rcsuhing system of equations can be 



6i = 1 - 62, 61 = 1 - 7 - ^2, 

6202 = 1/2, (1-7-62)01 +62C2 = 1/2-7, 

(1-62)01+6202 = 1/2, 6252 = 1/2-7- 

In order to solve system (|8.1ip we can compute the coefficients as follows 

62 = 1/(252), 62 = (1 - 27)/(252), 

and 

62(02 - Cl) = 1/2 - 7 - 01 + C17, 
62(02-01) = 1/2 -Cl. 

Substituting 62 and 62 in (j8.12p we get 

(c2 - Cl) _ 1/2 - 7 - Cl + ci7 



(8.11) 



(8.12) 



252 

(C2 - Cl) 

252 



1-27 
1/2 -01. 



Now, comparing and equating the two expressions we have 7 = and it is impossible because the matrix 
A is invertible. 



8.3. Second and third order IMEX schemes. 

1. Second order IMEX schemes: 
- ARS(2,2,2) scheme, [J 
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7 










7 





7 





1 


S 


1 - 


,5 





1 





1-7 


7 




S 


1 - 


,5 










1-7 


7 



with 7 = (2 - V2)/2 and (5 = 1 - 1/(27). 
- SSP2-(3,3,2) scheme, [36] 















1/4 


1/4 








1/2 


1/2 








1/4 





1/4 
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1/2 


1/2 
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1/3 


1/3 


1/3 




1/3 


1/3 


1/3 




1/3 


1/3 


1/3 



2. Third order IMEX schemes: 
- ARS(4,4,3) scheme, [I] 
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- BPR(3,5,3) scheme 
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